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ABSTRACT 

Electric  propulsion  thrusters  are  replacing  small  chemical  thrusters  used  for  spacecraft  control  and  orbital 
maneuvers.  These  thrusters  use  a  variety  of  mechanisms  to  convert  electrical  power  into  thrust  and  in  general 
provide  superior  specific  impulse  in  comparison  to  chemical  systems.  Electric  propulsion  has  been  under 
development  for  the  last  fifty  years,  and  almost  all  thrusters  are  designed  based  largely  on  experience  and 
experimentation.  The  present  article  considers  the  progress  made  in  numerical  simulation  of  electric 
propulsion  thrusters.  Due  to  the  wide  range  of  such  devices,  attention  is  restricted  to  electric  propulsion 
thruster  types  that  are  presently  in  use  by  orbiting  spacecraft.  The  physical  regimes  created  in  these  thrusters 
indicate  that  a  variety  of  numerical  methods  are  required  for  accurate  numerical  simulation  ranging  from 
continuum  formulations  to  kinetic  approaches.  Successes  of  numerical  simulation  models  are  demonstrated 
through  specific  examples.  It  is  concluded  that  numerical  simulations  can  be  expected  to  play  a  more 
prominent  role  in  the  design  and  evolution  of  future  electric  propulsion  thrusters. 


1.0  INTRODUCTION 


Electric  propulsion  technology  generates  thrust  primarily  from  electrical  energy  through  a  number  of  different 
mechanisms.  In  general,  electric  thrusters  provide  superior  performance  for  spacecraft  propulsion  in 
comparison  to  chemical  thrusters  due  to  their  significantly  higher  specific  impulses.  The  advantage  of  high 
specific  impulse  propulsion  can  be  understood  from  the  rocket  equation: 
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(1.1) 


where  M.  and  Mf  are  the  initial  and  final  spacecraft  masses,  respectively,  AUis  the  required  change  in 
spacecraft  velocity  for  the  mission,  g  is  the  gravitational  acceleration  (at  the  Earth’s  surface),  and  I  is  the 

specific  impulse.  Equation  (1.1)  indicates  that  to  maximize  useful  payload,  we  want  to  maximize  the  specific 
impulse.  While  propulsion  systems  based  on  chemical  propellants  are  limited  to  Isp  in  the  range  of  200-400 

sec,  electric  propulsion  systems  generate  Isp  in  the  range  of  several  hundred  to  several  thousand  seconds. 

Table  1  lists  typical  values  of  specific  impulse  and  thrust  associated  with  several  types  of  electric  propulsion 
systems.  In  addition  to  superior  propellant  mass  efficiency,  electric  propulsion  is  ideal  for  providing  very  low 
impulse  bits  required  by  high-resolution  observing  missions  and  by  spacecraft  constellations.  There  are  two 
categories  of  activities  involving  numerical  modelling  of  electric  thrusters.  The  first  concerns  modelling  of 
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the  thruster  device  and  its  acceleration  processes.  The  primary  goal  of  these  activities  is  to  aid  in  the 
fundamental  understanding  of  the  physical  mechanisms  of  the  thruster  with  the  objective  of  improving  its 
performance  in  terms  of  propulsion  efficiency  and  operational  lifetime.  The  second  area  of  modelling  activity 
concerns  the  plumes  produced  by  electric  thrusters.  Detailed  information  on  the  plumes  is  required  for  safe 
integration  of  the  thruster  onto  a  spacecraft.  Modelling  plays  an  important  role  here  because  it  is  problematic 
and  expensive  to  reproduce  the  in-orbit  space  environment  using  ground-based  laboratory  facilities.  Device 
modelling  also  plays  an  important  role  in  plume  simulations  by  providing  accurate  boundary  conditions  at  the 
thruster  exit.  The  present  article  is  focused  on  the  numerical  simulation  of  the  flows  in  electric  propulsion 
thrusters. 

While  there  is  ample  motivation  for  the  use  of  numerical  modelling  in  electric  propulsion,  much  of  the 
thruster  design  and  spacecraft  integration  work  continues  to  be  experimental  and  empirical.  This  is  largely 
due  to  the  technical  challenges  of  creating  physically  accurate  models  of  these  complex  systems.  In  this 
article,  the  current  state  of  the  art  for  numerical  modelling  of  the  device  flows  of  a  variety  of  electric  thrusters 
is  reviewed.  Due  to  space  limitations,  the  focus  is  restricted  to  thruster  types  that  are  presently  operating  in 
space.  Thus,  the  thrusters  considered  are  resisto-jets,  arc -jets,  gridded  electro-static  ion  thrusters,  and  Hall 
effect  thrusters.  Several  other  thruster  types  are  therefore  omitted  including  magnetoplasmadynamic  (MPD) 
thrusters,  pulsed  plasma  thrusters,  thrusters  based  on  electro-magnetic  wave  sources,  field  emission  electric 
propulsion  (FEEP),  colloid  thrusters,  and  laser-based  thrusters.  The  fundamental  physical  characteristics  of 
the  thrusters  considered  are  discussed  to  illustrate  the  types  of  numerical  approaches  required  to  accurately 
model  the  most  important  phenomena.  It  will  become  clear  that  some  thrusters  can  be  modelled  using  a 
continuum  formulation  while  others  require  a  kinetic  approach.  Wherever  appropriate,  the  strengths  and 
weaknesses  of  the  existing  models  for  these  thrusters  are  discussed. 

Finally,  the  future  prospects  for  improvement  in  numerical  modelling  of  electric  thrusters  are  considered.  As 
with  many  technological  areas,  the  potential  for  numerical  modelling  in  the  design,  optimization,  and 
integration  of  electric  thrusters  is  significant.  However,  the  physical  complexity  of  many  of  these  thrusters 
makes  the  goal  of  developing  accurate  numerical  models  one  that  still  requires  much  fundamental  research. 
Specific  areas  where  advances  are  required  in  the  modelling  are  discussed. 


2.0  NUMERICAL  METHODS 

The  focus  of  this  article  is  on  numerical  methods  used  to  model  the  flow  of  gas  and  plasma  through  electric 
propulsion  devices.  Discussion  of  the  numerical  analysis  of  other  aspects  of  thruster  design,  such  as  thermal 
and  structural  processes,  is  omitted  here.  There  are  two  fundamental  types  of  numerical  approaches  for 
modelling  these  gas  and  plasma  flows:  (1)  continuum  or  fluid  methods,  and  (2)  kinetic  methods.  The  choice 
of  method  for  modelling  a  specific  thruster  should  be  dictated  by  the  physical  characteristics  of  the  flow  in  the 
device,  and  by  the  level  of  accuracy  required  from  the  simulation.  As  we  shall  see,  there  is  a  wide  range  of 
conditions  in  different  types  of  thrusters.  In  general  this  means  that  different  methods  and  computer  codes 
must  be  developed  for  each  of  the  different  thrusters.  Before  considering  models  for  the  various  thrusters,  let 
us  first  review  the  most  important  numerical  methods  that  are  employed. 

2.1  Continuum/Fluid  Methods 

The  gas  or  plasma  may  be  considered  as  a  continuum  or  fluid  when  length  scales  associated  with  the  most 
important  physical  phenomena  (spatial  gradients  of  flow  properties,  distance  between  collisions,  charge 
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separation  distance)  are  small  in  comparison  to  the  size  of  the  thruster.  In  general,  these  conditions  are 
encountered  in  the  relatively  high  density  devices  such  as  resisto-jets  and  arcjets.  Continuum  methods  are 
usually  based  on  a  set  of  partial  differential  equations  describing  conservation  of  mass  (or  charge), 
momentum,  and  energy.  For  example,  the  viscous  flow  encountered  in  a  resisto-jet  may  be  modelled  using 
the  following  equation  set  (often  referred  to  as  the  Navier-Stokes  equations): 


dp  dpuj 

Continuity:  —  + - -  =  Sc 

dt  3c  j 


Momentum: 


dpp  dpuiiti 
dt  ckj 


^  dKk 


+  pF{+  SM. 


Energy: 


dpe  dpu  e  dp  d  r 
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(2.2) 


(2.3) 


where  p  is  the  mass  density,  uL  is  a  velocity  component,/?  is  the  pressure,  e  is  the  total  specific  internal  energy, 
Ty  is  the  shear  stress  tensor,  qk  is  the  heat  flux  vector,  is  an  external  body  force,  and  Sc,  SM  ,  and  SE  are 

source  terms  caused  by  chemical  reactions  in  the  mass,  momentum,  and  energy  equations,  respectively.  First- 
order  models  for  the  shear  stress  and  heat  flux  vector  use  temperature-dependent  transport  coefficients 
(viscosity  and  thermal  conductivity).  For  plasma  flows,  additional  equations  supplement  the  above  set.  These 
additional  equations  are  highly  specific  to  the  particular  thruster  type.  In  general,  it  is  necessary  to  include  an 
additional  set  of  conservation  equations  for  the  electrons  and  to  add  terms  to  the  heavy  particle  momentum 
equation  to  account  for  the  electro-magnetic  fields  and  friction  with  the  electrons.  Equation  sets  of  this  type 
can  be  solved  using  a  wide  variety  of  numerical  approaches.  In  finite  difference  formulations,  each  partial 
derivative  is  evaluated  to  some  order  using  Taylor  expansions.  In  finite  volume  formulations,  the  integral 
forms  of  the  conservation  equations  are  solved.  The  issues  regarding  these  formulations  are  discussed  in  a 
variety  of  texts,  e.g.  Tannehill,  Anderson  &  Pletcher  [1],  Fletcher  [2]. 


2.2  Kinetic  Methods 

The  flows  in  many  electric  propulsion  systems  fall  into  the  rarefied  regime  in  which  collision  and  plasma 
length  scales  are  similar  to  or  even  larger  than  the  size  of  the  thruster.  In  such  cases,  the  flow  is  not  well 
represented  by  a  continuum  formulation,  and  instead,  a  molecular,  kinetic  approach  must  be  undertaken.  Let 
us  discuss  separately  the  numerical  simulation  of  rarefied  gas-dynamic  and  plasma-dynamic  effects. 


2.2.1  Gas  Dynamics 

The  Knudsen  number  is  defined  as  the  ratio  of  the  mean  free  path  (the  average  distance  travelled  by  each 
molecule  between  successive  collisions)  and  the  characteristic  length  of  the  flow  (e.g.  the  diameter  of  a 
thruster).  In  general,  the  rarefied  flow  regime  is  defined  to  exist  for  Knudsen  numbers  above  0.01.  The 
continuum  approach  fails  under  rarefied  conditions  because  there  is  an  insufficient  rate  of  collisions  to 
maintain  the  velocity  distribution  functions  anywhere  near  to  the  small  departures  from  the  Maxwellian 
equilibrium  form  implicitly  assumed  in  continuum  formulations.  While,  in  principle,  the  Boltzmann  equation 
of  dilute  gas  dynamics  should  be  solved  under  high  Knudsen  number  conditions,  this  turns  out  to  be  a 
formidable  mathematical  and  computational  task.  The  direct  simulation  Monte  Carlo  method  (DSMC)  of  Bird 
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[3]  is  a  well-developed,  highly  successful  numerical  technique  for  simulating  rarefied  gas  flows.  In  this 
technique,  no  attempt  is  made  to  solve  any  specific  governing  equation.  Instead,  a  direct  simulation  of  the 
molecular-level  gas  dynamics  is  performed.  In  the  DSMC  technique,  a  large  number  of  model  particles  is 
simulated.  Each  model  particle  represents  a  much  larger  number  of  real  atoms  or  molecules.  The  model 
particles  possess  unique  velocity  components.  A  key  aspect  of  the  DSMC  technique  is  that  particle  motion  is 
decoupled  from  particle  collisions.  This  is  only  reasonable  if  the  time  step  employed  in  the  simulation  is 
significantly  smaller  than  the  mean  time  between  collisions.  Thus,  during  each  iteration  of  the  DSMC 
algorithm,  all  the  particles  are  first  moved  according  to  the  product  of  the  time  step  and  their  individual 
velocity  vectors.  Any  interactions  between  the  particles  and  boundaries  are  then  processed.  The  particles  are 
next  collected  into  cells  for  computation  of  collisions.  Within  each  cell,  the  positions  of  the  particles  are 
ignored.  Instead,  the  particles  are  paired  up  at  random  and  collision  probabilities  are  computed  for  each  pair 
to  determine  whether  a  collision  occurs.  This  statistical  approach  is  only  accurate  when  the  size  of  the 
computational  cell  is  of  the  order  of  a  mean  free  path.  The  collision  probability  is  proportional  to  the  product 
of  the  collision  cross  section  and  the  relative  velocity  of  the  colliding  particles.  Many  different  types  of 
collision  phenomena  can  be  simulated  using  the  DSMC  technique  including  momentum  exchange,  charge 
exchange,  internal  energy  exchange,  and  chemical  reactions.  Average  flow  properties  such  as  density, 
velocity,  and  pressure  are  obtained  by  time-averaging  of  the  particle  properties.  Unsteady  flows  may  be 
simulated  by  averaging  over  small  periods  of  time,  or  by  ensemble  averaging.  A  related  method,  the  Monte 
Carlo  Collision  (MCC)  technique  [4],  is  also  often  used  in  electric  propulsion  simulations.  This  approach  is 
useful  in  hybrid  particle-fluid  simulations  to  determine  the  effects  of  collisions  of  electrons  (modelled  as  a 
fluid)  on  the  heavy  species  that  are  represented  as  particles. 

2.2.2  Plasma  Dynamics 

Similar  to  gas  flow,  a  plasma  cannot  be  described  accurately  by  a  continuum  formulation  under  high  Knudsen 
number  conditions.  However,  continuum  models  can  still  yield  useful  solutions  in  cases  where  the  underlying 
collision  effects  are  dominated  by  the  plasmadynamics.  For  electric  propulsion  devices,  the  most  important 
behaviour  to  simulate  is  the  plasmadynamic  effects  on  the  momentum  of  the  heavy  ions  since  these  are  the 
main  source  of  thrust.  The  Particle  In  Cell  (PIC)  method,  see  Birdsall  and  Langdon  [5],  is  a  particle  technique 
for  simulating  plasmas.  In  most  electric  propulsion  systems,  ions  are  treated  as  particles  using  PIC  and  the 
electrons  are  treated  using  a  continuum,  fluid  approach.  This  approach  is  usually  justified  as  the  time  scales 
associated  with  the  much  lighter  electrons  are  orders  of  magnitude  smaller  than  those  for  the  heavy  ions.  In 
addition,  even  if  the  velocity  distribution  function  of  the  electrons  is  not  Maxwellian,  it  is  usually  compact 
allowing  meaningful  interpretation  of  continuum  quantities  such  as  temperature  and  pressure.  Plasma  physics 
has  been  described  as  the  art  of  approximation,  and  this  is  an  important  aspect  of  applying  the  PIC  method  to 
electric  propulsion  flows.  The  underlying  philosophy  should  always  be  to  take  the  simplest  approach  that 
provides  adequate  results.  Wherever  possible,  the  assumption  of  charge  quasi-neutrality  is  employed.  This 
approach  allows  the  spatial  distribution  of  the  ions  to  give  the  spatial  distribution  of  the  electrons.  This  is  a 
good  assumption  provided  the  Debye  length,  the  distance  over  which  electrons  and  ions  are  shielded  from  one 
another,  is  much  smaller  than  the  length  scale  of  the  thruster.  For  a  typical  thruster  electron  temperature  of 
several  electron-volts,  the  assumption  of  quasi-neutrality  is  good  down  to  plasma  densities  of  about  1014  m"3 
that  is  much  lower  than  values  encountered  in  most  electric  propulsion  thrusters  (see  values  listed  in  Table  1) 
except  for  the  inter- grid  regions  of  ion  thrusters.  A  variety  of  approximations  of  the  electron  momentum  and 
energy  equations  are  invoked  for  different  types  of  electric  propulsion  thrusters.  These  will  be  discussed  for 
each  thruster  type.  The  electric  field  is  usually  obtained  from  spatial  differentiation  of  the  plasma  potential 
that  in  turn  is  derived  from  the  electron  momentum  equation.  This  field  is  used  to  change  the  velocity  vector 
of  each  ion  particle  in  the  PIC  simulation  according  to  Newton's  2nd  Faw: 
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at 

where  is  the  ion  mass,  v  is  the  ion  velocity  vector,  qt  is  the  ion  charge,  and  E  is  the  electric  field.  The 
steps  in  the  PIC  algorithm  are  to  obtain  the  plasma  density  from  the  spatial  distribution  of  the  ions,  calculate 
the  electric  field,  accelerate  the  ions,  and  finally  move  the  ions.  Note  that  the  DSMC  calculation  of  ion 
collisions  is  very  easily  introduced  into  the  PIC  algorithm  after  the  movement  step. 


3.0  THRUSTER  MODELS 

In  the  following  sections,  a  summary  is  provided  of  the  past  history  and  current  status  of  numerical  modelling 
of  four  different  types  of  electric  propulsion  thrusters  presently  operating  in  space. 

3.1  Resisto-Jets 

A  resisto-jet  is  classified  as  an  electrothermal  device.  It  operates  by  flowing  gaseous  propellant  over  an 
electrical  element  that  raises  the  gas  temperature.  The  heated  gas  is  accelerated  gas  dynamically  through  a 
converging-diverging  nozzle  to  produce  thrust.  As  of  May  2010,  there  were  more  than  100  operational 
spacecraft  using  hydrazine  resisto-jets  in  space.  Resisto-jets  typically  have  good  efficiencies  of  70-80%.  The 
key  performance  loss  mechanisms  concern  the  development  of  relatively  large  viscous  boundary  layers  along 
the  nozzle  wall  (caused  by  the  relatively  low  gas  density),  and  frozen  flow  losses  caused  by  the  freezing  of 
energy  into  internal  modes  such  as  rotation  and  vibration. 

The  physical  characteristics  of  a  typical  resisto-jet  involve  a  number  density  at  the  nozzle  exit  of  1021  m'3,  a 
velocity  of  1-5  km/s,  and  an  exit  diameter  of  a  few  centimetres.  These  values  give  a  Knudsen  number  at  the 
nozzle  exit  of  about  0.1  that  indicates  the  gas  flow  involves  transition  from  continuum  flow  at  the  nozzle 
throat,  to  rarefied  flow  at  the  nozzle  exit.  Numerical  simulations  of  resisto-jets  have  employed  both 
continuum  and  kinetic  methods.  For  example,  Kim  [6]  solved  the  Navier-Stokes  equations  using  a  steady, 
finite-volume  formulation.  Full  DSMC  computations  of  the  nozzle  flows  of  cold  and  hot  resisto-jet  flows  were 
performed  by  Boyd  et  al.  [7,8].  Figure  1  shows  Mach  number  contours  computed  using  DSMC  in  the 
diverging  nozzle  of  a  cold-flow  hydrogen  resisto-jet  [8].  Such  computations  typically  are  begun  at  the  nozzle 
throat  assuming  sonic  conditions.  The  contours  clearly  illustrate  the  thick,  viscous  boundary  layer  formed 
along  the  nozzle  wall.  Figure  2  provides  a  comparison  between  experimental  measurements  (obtained  using 
Coherent  Anti  Raman  Scattering,  CARS)  and  DSMC  simulation  for  the  velocity  along  the  axis  of  the 
hydrogen  resisto-jet  [8].  Simulations  and  experimental  measurements  are  again  compared  in  Fig.  3  for 
specific  impulse.  The  simulations  agree  within  5%  of  the  measured  data.  A  comparison  between  data 
obtained  with  the  Navier-Stokes  equations,  with  the  DSMC  method,  and  with  experimental  measurements  of 
pitot  pressure  was  conducted  by  Boyd  et  al.  [7].  It  was  found  that  the  DSMC  computations  offered  superior 
agreement  with  the  measurements  than  that  obtained  with  the  CFD  analysis. 

As  demonstrated  by  the  results  shown  in  Figs.  2,  3,  the  accuracy  of  numerical  simulations  for  resisto-jet  flows 
is  high  due  to  the  relative  simplicity  of  the  physical  phenomena  that  occur  in  this  device.  Thus,  computational 
analysis  can  be  used  with  confidence  to  aid  in  the  design  of  new  resisto-jets.  This  is  the  approach  taken  for  the 
Free-Molecule  Microresistojet  developed  for  low  thrust  levels  that  has  been  designed  extensively  by 
Ketsdever  et  al.  [9]  using  the  DSMC  technique.  Computed  values  for  performance  parameters  were  found  to 
agree  with  experimental  measurements  to  within  2-4%  for  different  propellants. 
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Although  many  resisto-jets  are  presently  operating  in  space,  interest  in  further  development  of  this  thruster 
type  has  declined  since  significantly  superior  performance  are  offered  by  the  other  forms  of  electric 
propulsion  thrusters.  Specialized  resisto-jet  systems  have  been  developed  and  tested  in  recent  years  that  use 
water,  xenon,  and  ammonia  as  propellants  for  small  spacecraft  [10,1 1,12]. 

3.2  Arc  jets 

An  arcjet  thruster  is  also  an  electrothermal  device.  It  operates  by  flowing  gas  through  an  electrical  arc  formed 
between  a  cathode  placed  near  the  throat  of  a  nozzle,  and  the  nozzle  itself  that  acts  as  the  anode.  The  arc 
transfers  energy  to  the  gas  through  the  Ohmic  (Joule)  heating  mechanism,  and  the  hot  gas  is  subsequently 
accelerated  through  a  diverging  nozzle.  As  of  May  2010,  there  were  more  than  40  operational  spacecraft 
using  hydrazine  arcjets  in  space.  The  efficiency  of  an  arcjet  is  only  about  30%  and  is  significantly  lower  than 
for  a  resisto-jet  due  to  more  significant  frozen  flow  losses.  These  losses  result  from  a  variety  of 
nonequilibrium  collisional  phenomena  that  occur  in  arcjet  nozzle  flows  including  rotational  and  vibrational 
energy  relaxation,  dissociation  of  molecules,  and  ionization  of  atoms.  As  illustrated  in  Table  1,  however, 
arcjets  do  offer  significantly  higher  specific  impulse  than  resisto-jets. 

Arcjets  developed  for  space  propulsion  typically  operate  at  power  levels  of  1  to  10  kW  using  hydrazine 
propellant,  and  have  physical  characteristics  that  are  similar  in  density  to  resisto-jets,  although  the  velocity  is 
usually  higher  resulting  in  the  higher  specific  impulses.  Most  arcjet  device  models  are  based  on  a  continuum 
approach  in  which  the  continuum  conservation  equations  are  solved  together  with  appropriate  transport 
coefficients.  Butler  and  King  [13]  developed  a  two-fluid  model  of  a  10  kW  hydrogen  arcjet  that  used  separate 
continuity,  momentum,  and  energy  equations  for  the  heavy  particles  and  the  electrons.  The  fluid  equations 
were  augmented  with  basic  plasma  relations  such  as  Faraday's  law,  Maxwell's  equation  for  the  magnetic  field, 
and  Ohm's  Law  for  Joule  heating.  In  a  similar  approach,  Miller  and  Martinez-Sanchez  [14]  solved  a 
continuity  equation  for  each  species,  momentum  equations  for  the  axial,  radial,  and  azimuthal  coordinates, 
and  energy  equations  for  the  electrons  and  the  heavy  particles  that  included  Ohmic  heating,  and  losses  due  to 
electron  collisions  and  radiation.  In  addition,  dissociation  and  ionization  reactions  for  the  hydrogen  propellant 
were  included.  These  models  typically  offer  agreement  within  10%  of  measured  performance  data.  A  similar 
arcjet  model  was  presented  by  Megli  et  al.  [15]  for  hydrazine  propellant  by  implementing  a  finite  rate, 
nonequilibrium,  chemical  kinetics  mechanism  for  the  N2-H2  system.  Models  developed  for  high  power  arcjets 
(100  kW)  by  Auweter-Kurtz  et  al.  [16]  add  magnetic  effects  to  the  momentum  and  energy  equations,  and  a 
separate  energy  equation  for  the  vibrational  energy  of  molecules.  One  example  of  a  particle-based  kinetic 
description  of  arcjets  was  developed  by  Boyd  [17,18]  in  which  the  usual  DSMC  technique  for  a  multi-species, 
reacting  gas,  is  augmented  by  a  kinetic  model  for  Ohmic  heating.  Continuum  and  DSMC  methods  were 
applied  to  investigate  the  effects  of  nozzle  geometry  on  arcjet  design  by  Hatta  and  Aso  [19]. 

The  accuracy  of  numerical  simulations  for  arcjet  flows  is  much  less  than  that  for  resisto-jets  due  to  the 
significantly  wider  range  of  physical  processes  that  are  involved.  The  reduction  in  accuracy  stems  mainly 
from  a  lack  of  information  on  some  of  these  physical  processes  such  as  molecular  transport  coefficients  in  a 
nonequilibrium  plasma.  Figure  4  (taken  from  Boyd  [18])  shows  the  decoupled  flow  domain  in  which  a 
continuum  approach  is  employed  up  to  the  nozzle  throat  that  is  used  as  a  start  line  for  the  DSMC  analysis. 
The  nozzle  wall  temperatures  are  also  taken  from  the  continuum  model.  Figures  5  and  6  show  comparisons  of 
axial  and  radial  velocity  profiles  in  a  1.4  kW  hydrogen  arcjet.  The  measured  data  were  obtained  using  Laser 
Induced  Fluorescence.  The  agreement  between  the  predicted  and  measured  data  is  quite  good.  Comparison 
between  the  DSMC  simulation  (with  and  without  Ohmic  heating  Q)  and  measurement  of  specific  impulse  is 
shown  in  Fig.  7  where  the  agreement  between  the  datasets  is  within  10%. 
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Similar  to  resisto-jets,  although  many  arcjets  are  presently  operating  in  space  on  communications  satellites, 
interest  in  further  development  of  this  thruster  type  has  declined  since  significantly  better  performance  are 
offered  by  ion  and  Hall  thrusters.  Notable  exceptions  include  development  of  arcjets  for  lunar  and  planetary 
exploration  [20,  21].  In  addition,  similar  to  other  EP  devices,  arcjet  designs  have  been  scaled  down  for  small 
satellite  applications  [22]. 

3.3  Ion  Thrusters 

A  gridded  ion  thruster  is  an  electrostatic  device  that  operates  by  flowing  gas  propellant  (usually  xenon  or 
krypton)  into  a  discharge  chamber  where  a  hollow  cathode  provides  electrons  that  ionize  the  gas.  The 
resulting  plasma  then  drifts  towards  a  set  of  grids  (usually  two:  a  screen  grid  and  an  acceleration  grid)  across 
which  a  large  potential  difference  is  applied  in  order  to  accelerate  the  ions  to  very  high  velocity.  The  ions 
flow  through  the  grids  via  a  large  number  of  small  apertures  each  of  which  focuses  the  ions  and  is  therefore 
called  the  ion  optics.  A  schematic  diagram  of  an  ion  thruster  is  provided  in  Fig.  8  that  shows  the  main 
components.  In  the  following,  the  status  of  modelling  is  reviewed  of  the  major  ion  thruster  components, 
namely  the  discharge  chamber ,  the  hollow  cathode ,  and  the  ion  optics.  As  of  May  2010,  there  were  more  than 
30  operational  spacecraft  using  xenon  gridded  ion  thrusters  in  space.  Notable  applications  of  ion  thrusters 
include  Deep  Space- 1  (NASA),  Dawn  (NASA),  Hayabusa  (JAXA)  and  the  XIPS  system  used  on  Boeing 
communications  satellites.  The  efficiency  of  an  ion  thruster  is  typically  about  70%.  The  primary  loss 
mechanisms  are  the  propellant  ionization  and  beam  divergence.  Another  key  issue  in  the  use  of  this  type  of 
electric  propulsion  system  is  thruster  lifetime.  Some  of  the  ions  accelerating  through  the  grids  can  impinge  on 
the  grids  at  energies  that  are  sufficiently  large  to  cause  sputtering  of  the  grid  material.  After  several  thousand 
hours  of  operation,  the  grids  can  become  completely  eroded  away  leading  to  failure  of  the  thruster. 

3.3.1  Discharge  Chamber 

A  key  step  in  ion  thruster  design  involves  achieving  uniform  plasma  conditions  in  the  discharge  chamber  and 
this  is  largely  controlled  by  employing  a  magnetic  field  topology  that  constrains  the  electron  trajectories  so  as 
to  enhance  ionization.  A  complete  model  of  the  discharge  chamber  of  an  ion  thruster  was  first  developed  by 
Arakawa  and  Ishihara  [23]  using  a  combination  of  a  finite  element  fluid  approach  and  a  Monte  Carlo  method 
for  the  primary  electrons  that  originate  from  the  cathode.  Ions  were  treated  as  a  fluid,  electro-static  fields  were 
ignored,  and  secondary  electrons  (those  created  in  the  ionization  process)  were  omitted.  In  a  more  advanced 
approach,  Jugroot  and  Harvey  [24]  applied  a  combination  of  the  DSMC  and  fluid  plasma  methods  to  analyze 
the  discharge  chamber  of  the  UK-T5  ion  thruster. 

Wirz  and  Katz  [25]  describe  the  development  of  a  comprehensive  axially  symmetric  discharge  chamber  model 
that  treats  the  primary  electrons  as  particles  using  PIC,  while  the  ions  and  secondary  electrons  are  described 
using  a  fluid  description.  The  neutral  gas  is  modelled  using  sources  and  sinks  based  on  expressions  from 
kinetic  theory.  An  even  more  comprehensive  modelling  approach  is  described  by  Mahalingam  and  Menart 
[26]  in  which  ions,  primary  electrons,  and  secondary  electrons  are  all  treated  as  particles  using  PIC  while  a 
uniform  neutral  gas  distribution  is  assumed.  Stueber  [27]  describes  the  development  of  a  three-dimensional, 
particle-based  simulation  technique  for  mapping  the  neutral  atom  distribution  that  is  to  be  connected  with  a 
three-dimensional,  fluid-based  description  of  the  electrons  [28]. 

Recent  progress  has  focused  on  detailed  and  numerically  efficient  modelling  of  electron  dynamics. 
Mahalingam  and  Menart  [29]  describe  application  of  their  numerical  tools  to  investigate  the  performance  of 
different  magnetic  field  topologies  in  confining  electrons.  The  same  authors  present  a  comprehensive 
assessment  of  their  computational  techniques  by  comparison  of  numerical  and  experimental  results  for  the 
NASA  NSTAR  ion  thruster  [30],  that  was  used  on  Deep  Space-1.  While  the  computed  ion  beam  current 
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agrees  with  experiments  within  1%,  the  computed  discharge  current  differs  from  the  measurements  by  more 
than  20%.  Recently,  a  fully  kinetic  approach  (PIC-MCC)  has  been  presented  by  Mahalingham  et  al.  [31]  for 
computation  of  the  electric  field  in  an  ion  thruster  discharge  chamber.  As  with  other  thruster  technologies,  the 
computational  models  developed  for  meso-scale  devices  are  also  being  applied  to  understand  the  plasma 
processes  in  smaller  scale  systems  [32]. 

3.3.2  Cathodes 

Two  hollow  cathodes  are  employed  in  a  typical  ion  thruster  system.  One  cathode  is  located  inside  the 
discharge  chamber  and  provides  the  source  of  ionizing  electrons.  The  second  cathode  is  located  external  to  the 
thruster  and  provides  a  source  of  electrons  to  neutralize  the  ion  beam  exiting  the  thruster.  In  general,  there  is 
interest  in  modelling  these  cathodes  as  their  erosion  is  an  important  life-limiting  mechanism  for  ion  thrusters. 

Due  to  the  relatively  high  number  density  (1022  m'3)  and  small  physical  size  (1  cm)  of  a  hollow  cathode,  much 
of  the  internal  flow  lies  in  the  continuum  regime.  Rapid  plume  expansion,  however,  soon  takes  the  flow  into 
the  kinetic  regime.  Relatively  little  research  has  been  conducted  on  the  development  of  numerical  models  of 
cathodes.  Application  of  a  conventional  computational  fluid  dynamics  model  of  the  flow  inside  a  hollow 
cathode  was  described  by  Murray  et  al.  [33]  that  omitted  much  of  the  important  plasma  physics.  A  detailed 
one-dimensional  model  of  the  plasma  flow  inside  a  hollow  cathode  was  described  by  Katz  et  al.  [34]. 
Significant  efforts  have  been  made  by  Mikellides  et  al.  [35,  36]  to  develop  a  detailed  model  of  the  plasma 
flow  inside  the  hollow  cathode  and  in  the  near-field  plume  expansion  region  near  the  keeper.  A  multi-fluid 
continuum  approach  is  employed  and  assessed  in  detail  through  comparisons  with  experimental  measurements 
for  the  NSTAR  thruster  [37].  An  important  conclusion  from  this  work  is  a  requirement  for  further  theoretical 
investigations  on  the  role  of  plasma  turbulence  in  affecting  cathode  erosion.  Another  concern  with  cathodes  is 
the  depletion  of  the  insert  material  as  this  reduces  the  ability  of  the  device  to  generate  electrons.  An 
axisymmetric  numerical  model  of  the  depletion  process  in  a  barium  oxide  cathode  is  described  by  Coletti  and 
Gabriel  [38]  and  assessed  through  direct  comparisons  with  measurements  obtained  for  several  different 
cathodes  [39]. 

3.3.3  Ion  Optics 

A  large  number  of  computer  models  for  ion  thruster  ion  optics  performance  have  been  developed  over  the 
years.  Most  of  these  models  rely  on  the  Particle  in  Cell  (PIC)  approach  for  simulating  the  rarefied  plasma 
conditions.  In  addition,  almost  all  of  these  codes  consider  the  plasma  flow  through  a  single  aperture  of  the 
grids  due  to  the  geometric  complexity  associated  with  a  real  thruster;  for  example,  the  NSTAR  thruster  has 
about  15,000  apertures. 

Arakawa  and  Ishihara  [23]  combined  their  discharge  model  with  a  grid  ion  optics  model  based  on  Poisson's 
equation,  the  continuity  equation,  and  Newton's  2nd  law.  Instead  of  using  the  PIC  technique,  a  flux  tube 
approach  was  adopted.  A  key  omission  from  that  model  is  the  effect  of  charge  exchange  collisions.  These 
interactions  involve  an  electron  being  transferred  from  a  slow  neutral  atom  to  a  fast  ion.  The  resulting  low 
energy  ions  are  very  responsive  to  the  electrostatic  fields  in  the  intergrid  region  and  it  is  these  ions  that  are 
most  likely  to  impinge  on  the  grids  causing  sputtering  and  erosion.  Peng  et  al.  [40]  describe  a  detailed 
simulation  scheme  employing  the  PIC  technique  and  a  Monte  Carlo  approach  for  simulating  the  charge 
exchange  events.  Further,  the  energy  distribution  of  ions  impinging  on  the  grids  was  collected  during  the 
simulation  and  used  to  predict  the  erosion  depths  on  the  grid  surface  for  1000  hours  of  thruster  operation. 
Bond  and  Latham  [41]  report  on  the  use  of  a  PIC  code  to  modify  the  optics  of  the  UK- 10  ion  thruster,  that 
employs  a  three-grid  configuration  in  which  a  third,  deceleration  grid  is  added.  Good  agreement  was  obtained 
between  experimental  measurements  and  the  code  predictions  for  the  total  current  collected  on  the 
acceleration  and  deceleration  grids  at  a  number  of  different  thruster  operating  conditions.  The  code  was 
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subsequently  employed  to  identify  a  modified  grid  geometry  that  would  provide  extended  thruster  lifetime. 
This  was  achieved  through  a  combination  of  a  smaller  diameter  for  the  acceleration  grid  holes  and  a  thicker 
acceleration  grid.  In  this  way,  the  predicted  lifetime  was  increased  by  almost  a  factor  of  two. 

Later  numerical  simulation  studies  of  ion  optics  continued  to  refine  and  develop  these  techniques  in  different 
ways.  For  example,  Boyd  and  Crofton  [42]  introduced  multiply  charged  ions  and  modelled  collisions  in  detail 
using  the  DSMC  technique.  A  typical  computational  domain  is  shown  in  Fig.  9.  The  plasma  from  the 
discharge  chamber  enters  the  flow  domain  along  the  left  hand  edge,  labeled  (1).  The  ions  flow  through  the 
three-grid  system  of  this  thruster,  biased  at  the  voltages  indicated,  and  exit  the  flow  domain  in  the  downstream 
plume  region,  labeled  (2).  In  their  study,  Boyd  and  Crofton  [42]  considered  several  thruster  operating 
conditions  of  the  UK- 10  thruster  and  excellent  agreement  with  measurements  was  obtained  for  acceleration 
grid  current  as  shown  in  Fig.  10.  Nakayama  and  Wilbur  [43]  describe  a  simplified  grid  optics  analysis 
approach  that  allows  thousands  of  parametric  studies  to  be  performed  in  just  a  few  hours.  Wang  et  al.  [44] 
described  a  three-dimensional  code  that  was  applied  to  model  erosion  on  the  NSTAR  thruster.  In  their  work, 
it  was  realized  that  relatively  long  computational  domains  downstream  of  the  acceleration  grid  are  required  in 
order  to  capture  all  of  the  back- flowing  charge  exchange  ions.  Figure  1 1  shows  a  comparison  of  the  model 
with  data  obtained  in  an  8,200  hour  life-test  of  the  thruster.  The  code  is  remarkably  successful  in  predicting 
the  pits-and-grooves  erosion  patterns  measured  experimentally  on  the  downstream  face  of  the  acceleration 
grid.  The  level  of  agreement  is  shown  quantitatively  in  Fig.  12  where  the  erosion  profiles  on  the  grid  are 
compared  along  a  line  joining  two  adjacent  apertures. 

Analysis  of  a  single  aperture  of  grid  optics  can  be  considered  as  one  of  the  most  mature  areas  of  numerical 
simulation  within  electric  propulsion.  However,  there  remains  the  important  question  of  how  to  integrate  the 
results  for  one  aperture  across  the  entire  grid  face.  This  is  not  a  straightforward  procedure  as  the  discharge 
chamber  usually  creates  some  variation  in  ion  current  density  across  the  screen  grid  while  the  profile  of 
neutral  atom  properties  tends  to  be  much  more  uniform.  These  issues  have  been  studied  in  detail  by  Emhoff 
and  Boyd  [45]  for  the  NEXT  ion  thruster,  in  use  on  NASA’s  DAWN  mission.  The  approach  followed  in  [45] 
was  to  perform  several  simulations  representing  conditions  corresponding  to  several  different  apertures  across 
the  thruster.  The  results  from  these  simulations  were  then  integrated  across  the  thruster  face  to  obtain  overall 
performance  data.  In  the  case  of  the  NEXT  thruster,  this  approach  yielded  good  agreement  within  2%  of 
measurements  of  thrust,  specific  impulse,  and  beam  current  [45].  In  general,  the  accuracy  of  simulations 
associated  with  ion  thruster  optics  is  high  in  terms  of  propulsion  performance,  which  explains  why  modelling 
now  plays  a  key  role  in  the  development  of  this  technology.  Accurate  simulation  of  lifetime  related  properties 
such  as  acceleration  grid  current  and  erosion  is  more  problematic  and  requires  continued  research. 

3.4  Hall  Thrusters 

A  Hall  thruster  is  essentially  a  grid-less  electrostatic  acceleration  device  that  operates  by  flowing  gas 
propellant  (usually  xenon  or  krypton)  through  holes  in  an  anode  into  an  annular  acceleration  channel. 
Electrons  are  supplied  by  an  external  hollow  cathode.  An  applied  magnetic  field  with  a  primarily  radial 
component  is  used  to  partially  trap  the  electrons  in  a  desired  region  of  the  acceleration  channel  leading  to  a 
concentrated  ionization  zone.  There  are  two  main  configurations  of  Hall  thrusters.  The  so-called  stationary 
plasma  thruster  (SPT)  employs  a  relatively  long  acceleration  channel  (of  a  few  centimeters)  and  ceramic  wall 
insulator  materials  (e.g.  boron  nitride  or  silicon  carbide).  The  anode  layer  thruster  (TAL)  employs  a  shorter 
acceleration  channel  (a  few  millimeters)  and  metallic  wall  materials  (e.g.  molybdenum).  A  schematic  diagram 
of  an  SPT  Hall  thruster  is  provided  in  Fig.  13  that  shows  the  main  components.  Typically,  for  both  SPT's  and 
TAL's,  a  potential  difference  between  anode  and  cathode  of  several  hundred  volts  is  applied  to  accelerate  the 
ions  out  of  the  annular  channel.  Specific  impulses  in  the  range  of  1500  to  2500  seconds  are  achieved  at  an 


RTO-EN-AVT-162 


PAPER  NBR  -  9 


UNCLASSIFIED 


UNCLASSIFIED 


Simulation  of  Electric  Propulsion  Thrusters 


ORGANIZATION 


overall  efficiency  of  about  50-60%.  The  thrust  level  and  specific  impulse  make  Hall  thrusters  ideal  for 
station-keeping  tasks  on  communications  spacecraft.  As  of  May  2010,  there  were  about  18  operational 
spacecraft  using  Hall  thrusters  in  space.  Notable  applications  of  Hall  thrusters  include  the  SMART- 1  mission 
(ESA)  and  many  Russian  communications  satellites.  The  number  of  Hall  thrusters  operating  in  space  will 
increase  significantly  in  the  next  few  years  as  commercial  companies  in  the  US  and  Europe  have  already 
begun  implementing  Hall  thrusters  on  communications  satellites  for  station-keeping. 

There  are  several  aspects  of  Hall  thruster  performance  where  numerical  modelling  can  play  a  role  in 
identifying  improved  operating  conditions.  Two  of  the  key  performance  loss  mechanisms  involve  collisions 
of  ions  with  the  walls  of  the  acceleration  channel,  and  beam  divergence  (see  discussions  in  Kim  [46]).  Ions 
colliding  with  the  walls  lose  their  momentum  and  are  reflected  as  slow  neutral  atoms.  In  order  to  contribute  to 
thrust,  these  particles  must  be  ionized  and  accelerated  again,  and  this  clearly  represents  a  performance  loss.  In 
addition,  the  ions  can  impact  the  thruster  walls  at  energies  that  are  sufficiently  large  to  sputter  the  wall 
material,  and  this  represents  one  of  the  primary  lifetime-limiting  mechanisms  of  Hall  thrusters.  The 
divergence  angle  of  the  plume  of  a  Hall  thruster  is  significantly  larger  than  that  for  a  gridded  ion  thruster  that 
leads  to  greater  concerns  over  spacecraft  integration  as  well  as  greater  thrust  loss.  The  main  reasons  for  the 
increased  plume  divergence  are  the  topology  of  the  magnetic  field  that  leaks  out  into  the  near-field  plume 
region,  and  the  larger  electron  temperature  in  a  Hall  thruster  that  leads  to  increased  expansion  of  the  plume. 

The  SPT-100  Hall  thruster,  developed  in  the  former  Soviet  Union,  produces  a  plasma  density  of  about  1017  m'3 
and  a  velocity  of  about  16  km/s  in  an  acceleration  channel  with  an  annular  diameter  of  100  mm.  The  peak 
magnetic  field  strength  is  about  200  G  that  means  the  electrons  are  magnetized  and  the  ions  are  not.  These 
physical  characteristics  indicate  that  while  the  electrons  may  be  modelled  using  a  fluid,  continuum  approach, 
the  ions  and  neutrals  require  a  kinetic  formulation. 

Almost  all  numerical  modelling  studies  of  Hall  thrusters  have  considered  the  SPT  configuration.  One  of  the 
first  investigations,  by  Komurasaki  and  Arakawa  [47],  employed  a  steady,  two-dimensional  formulation  using 
a  fluid  electron  model  and  a  flux-tube  ion  model  (similar  to  that  used  in  ion  thruster  modeling).  This  approach 
was  used  to  investigate  basic  plasma  behaviour  and  thruster  performance  including  ion  loss  to  the  walls. 
There  are  several  important  limitations  in  the  approach  of  [47]  including  the  omission  of  the  neutrals  and  the 
assumption  of  steady  state  behaviour.  Experimental  measurements  of  Hall  thrusters  have  detected  a  variety  of 
oscillations  including  some  in  the  range  of  10-50  kHz  that  may  play  an  important  role  in  overall  device 
performance.  Such  oscillations  were  predicted  numerically  in  the  one-dimensional,  unsteady  study  of  Boeuf 
and  Garrigues  [48]  in  which  ions  were  modelled  using  the  kinetic  Vlasov  equation  and  the  neutrals  were 
included  using  a  constant  velocity  and  the  continuity  equation.  A  key  limitation  of  this  model  is  a 
phenomenological  treatment  of  wall  effects.  A  related  one-dimensional  model  in  which  the  ions  were 
modelled  using  a  PIC  approach  was  described  by  Garrigues  et  al.  [49],  in  which  good  qualitative  agreement 
was  obtained  with  experimental  data  for  SPT  performance  characteristics  for  both  xenon  and  krypton 
propellants,  as  shown  in  Fig.  14.  Several  alternative  one-dimensional  Hall  thruster  models  have  been 
successfully  developed  to  explore  fundamental  properties  of  Hall  thrusters.  For  example,  Fruchtman  et  al. 
[50]  used  a  fully  fluid  description  to  investigate  the  idea  of  controlling  the  electric  field  within  the  acceleration 
channel  using  absorbing  electrodes,  and  Ahedo  et  al.  [51]  employed  a  three-fluid  description  (electrons,  ions, 
and  neutrals)  to  investigate  effects  of  electron  pressure  and  back-flow  of  ions  to  the  anode.  Similar 
formulations  have  been  used  to  study  other  physical  phenomena  such  as  secondary  electron  emission  from  the 
walls  [52]  and  two-stage  Hall  thrusters  [53]. 

A  variety  of  two-dimensional,  unsteady  Hall  thruster  models  have  also  been  developed.  Fife  and  Martinez- 
Sanchez  [54]  used  a  self-consistent  PIC  model  for  ions  and  a  quasi-one-dimensional  fluid  model  for  electrons. 
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A  key  aspect  of  this  approach  is  to  sub-cycle  in  time  on  the  electron  properties  and  thus  the  model  is 
numerically  expensive.  This  same  model  was  used  by  Fife  et  al.  [55]  to  perform  a  detailed  analysis  of  the 
plasma-wall  interaction  phenomena  and  included  the  first  numerical  predictions  of  the  plasma  oscillations  also 
computed  by  Boeuf  and  Garrigues  [48].  A  similar  two-dimensional  model  was  reported  by  Koo  and  Boyd  [56] 
and  applied  to  the  AFRL/UM-P5  Hall  thruster.  Figure  15  illustrates  a  typical  computational  domain  for  these 
models  in  which  the  active  simulation  region  lies  between  a  virtual  anode  located  inside  the  thruster  and  a 
virtual  cathode  located  outside  the  thruster  in  the  near-field  plume  region.  In  a  sequence  of  investigations, 
Boeuf  et  al.  [57,58,59,60]  developed  and  extended  the  formulation  of  Fife,  and  used  their  model  to  investigate 
a  number  of  phenomena  including  the  role  of  electron  mobility  on  the  computational  results,  the  effects  of  the 
magnetic  field  on  thruster  performance,  and  analysis  of  a  two-stage  Hall  thruster.  The  models  reported  in  [54- 
60]  represent  the  state  of  the  art  in  modelling  Hall  thruster  performance.  The  accuracy  of  these  models  is 
found  to  depend  heavily  on  the  manner  in  which  electron  mobility  is  treated.  It  is  found  experimentally  that 
the  electron  current  obtained  in  the  operation  of  Hall  thrusters  is  higher  than  would  be  expected  from  classical 
modeling  of  electron  mobility  based  on  collisions  between  atoms  and  electrons.  Attempts  have  been  made  to 
explain  the  enhanced  mobility  through  electron  collisions  with  the  walls  of  the  thruster  (called  near-wall 
conductivity),  and  by  plasma  oscillations  or  turbulence  (called  Bohm  mobility).  Near-wall  conductivity  has 
been  analyzed  in  detail  by  Latocha  et  al.  [61].  Using  a  fully  kintic,  two-dimensional  approach,  Adam  et  al. 
[62]  show  that  plasma  turbulence  may  be  generated  under  the  high  ion-velocity  conditions  at  the  exit  of  a  Hall 
thruster  and  that  such  turbulence  can  lead  to  enhanced  electron  mobility.  Despite  these  efforts,  no  general 
mobility  model  has  yet  been  formulated  to  explain  all  the  phenomena  observed  in  thruster  operation  and  the 
issue  of  effective  electron  mobility  modelling  remains  a  critical  area  of  on-going  research  for  Hall  thruster 
simulations.  The  current  status  is  illustrated  by  Fig.  16  that  shows  profiles  of  plasma  potential  along  the  center 
of  the  annular  acceleration  channel  of  the  UM/AFRL  P5  Hall  thruster  [56].  It  is  found  that  the  experimental 
data  lies  between  the  two  simulations  that  employ  models  for  Bohm  mobility  and  wall-collision  mobility. 

The  hybrid  code  developed  by  Fife  and  Martinez-Sanchez  [54]  is  called  HPHall  and  this  code  has  served  as 
the  basis  for  continued  model  development.  The  most  recent  version  of  the  code,  called  HPHall2,  has  been 
advanced  significantly  by  Hofer  et  al.  [63,64,65].  In  Ref.  63,  new  models  for  propellant  injection  and  thruster 
wall  reflections  were  found  to  improve  agreement  of  computed  results  with  experimental  measurements.  In 
Ref.  64,  the  wall  sheath  model  of  Hobbs  and  Wesson  [66]  was  implemented  and  again  found  to  improve  the 
code  predictions.  Finally,  a  multiple-domain  electron  mobility  model  is  introduced  in  Ref.  65  and  found  to 
offer  the  modeller  increased  control  over  this  important  phenomenon. 

Efforts  are  continuing  to  develop  more  generalized  approaches  for  modelling  of  the  electrons  in  Hall  thrusters. 
A  shear-based  model  for  electron  transport  was  described  by  Scharfe  et  al.  [67]  but  found  to  give  agreement 
within  no  better  than  30%  with  experimental  measurements  of  thruster  performance.  A  fully  two-dimensional 
model  for  the  electron  dynamics  was  formulated  by  Escobar  and  Ahedo  [68].  Despite  the  lack  of  a  first 
principles  electron  mobility  model,  once  a  phenomenological  model  has  been  tuned  for  a  particular  type  of 
thruster,  numerical  simulations  can  play  a  role  in  helping  to  better  understand  the  complex,  coupled  plasma 
physics  mechanisms  present  in  Hall  thrusters.  For  example,  Scharfe  et  al.  [69]  use  measurements  of  various 
thruster  flow  quantities  to  assess  their  numerical  model  and  employ  it  to  study  effects  of  charge  exchange  and 
wall  collisions.  In  two  related  studies,  a  tuned  hybrid  Hall  thruster  model  was  employed  to  investigate  the 
potential  benefits  of  a  two-stage  device  in  which  the  ionization  and  acceleration  processes  are  physically 
separated  [70,71]. 

A  two-dimensional,  steady,  multi-fluid  formulation  has  been  developed  for  Hall  thrusters  by  Keidar  et  al.  [72]. 
These  simulations  require  run  times  of  minutes  rather  than  the  many  hours  required  by  the  unsteady,  hybrid 
fluid-particle  methods.  The  fully  fluid  methods  are  therefore  well-suited  to  the  rapid  analysis  of  critical 
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physical  phenomena  such  as  sheaths  [72]  and  magnetic  mirror  effects  [73].  The  rapid  simulation  times  offered 
by  fluid  formulations  also  make  them  very  well  suited  to  the  investigation  of  thruster  lifetime  issues.  For 
example,  the  change  in  Hall  thruster  wall  profiles  and  overall  thruster  performance  over  several  thousand 
hours  of  operation  of  the  SPT-100  Hall  thruster  were  investigated  by  Yim  et  al.  [74]  using  a  multi-fluid,  two- 
dimensional  model.  Profiles  of  the  predicted  time  evolution  of  the  eroded  walls  are  shown  in  Fig.  17  to  offer 
very  good  correspondence  to  the  experimentally  measured  behaviour.  Significantly  more  expensive  erosion 
studies  have  also  been  conducted  based  on  hybrid  simulations  of  the  discharge  plasma  [75,76]. 

A  two-dimensional,  unsteady,  multi-fluid  Hall  thruster  formulation  is  described  by  Mikellides  et  al.  [77]  that 
is  based  on  methods  employed  for  modelling  ion  thruster  hollow  cathodes  [36,37].  This  approach  offers 
several  advantages  over  the  hybrid-PIC  technique  including  a  fully  two-dimensional  electron  model,  the 
ability  to  extend  the  flow  domain  well  outside  the  thruster,  and  smooth  results  that  are  free  from  the  statistical 
fluctuations  inherent  in  particle  simulations.  The  method  has  been  applied  by  Mikellides  et  al.  [78]  to  explain 
unexpectedly  low  erosion  rates  measured  on  a  5  kW-class  Hall  thruster.  Of  course,  this  approach  still  requires 
input  of  an  electron  mobility  model  and  cannot  capture  all  of  the  nonequilibrium  effects  experienced  by  the 
discharge  ions. 

Anode  layer  Hall  thrusters  have  received  considerably  less  attention  in  terms  of  both  thruster  development  and 
numerical  simulation.  A  multi-fluid  model  of  a  high-power  bismuth  TAL  was  developed  by  Keidar  and  Boyd 
[79].  The  model  provided  good  agreement  with  measured  current-voltage  characteristics  and  was  later  further 
developed  for  application  to  a  two-stage  TAL  thruster  [80].  Due  to  the  relatively  small  size  of  the  acceleration 
channel  in  a  TAL  device,  it  can  be  argued  that  a  particle  approach  is  required  to  model  both  the  electrons  and 
the  ions. 

A  fully  kinetic,  two-dimensional  model  for  analysis  of  a  50-W  micro-TAL  Hall  thruster  was  formulated  by 
Szabo  and  Martinez-Sanchez  [81].  All  physical  phenomena  were  modeled  on  the  electron  time-scale.  Certain 
non-physical  assumptions  had  to  be  invoked  to  make  the  problem  numerically  tractable,  such  as  increasing  the 
electron-to-ion  mass  ratio,  and  increasing  the  permittivity  of  free  space.  Fully  kinetic  analyses  of  SPT-like 
thrusters  have  also  been  presented  by  Blateau  et  al.  [82]  and  Irishkov  et  al.  [83]  employing  similar  simplifying 
assumptions.  In  both  Refs.  82  and  83,  the  simulations  over  predict  the  measured  thrust  by  about  10%  and  the 
efficiency  by  about  50%.  Clearly,  there  is  still  significant  progress  required  to  improve  various  aspects  of 
fully  kinetic  simulation  of  Hall  thrusters.  Nevertheless,  this  is  a  likely  direction  for  future  research  as 
computer  power  continues  to  increase  allowing  more  accurate  simulation  of  the  electron  dynamics. 

Due  primarily  to  uncertainties  in  electron  mobility  modelling,  both  fluid  and  hybrid  simulations  schemes  are 
generally  only  capable  of  predicting  Hall  thruster  performance  (thrust  and  power)  within  about  10%. 
However,  once  a  simulation  method  has  been  tuned  using  a  phenomenological  mobility  model,  better 
accuracy  can  be  expected  over  a  limited  range  of  thruster  operating  conditions. 


4.0  SUMMARY 

The  successful  development  of  physically  accurate  numerical  methods  for  simulating  the  gas  and  plasma 
flows  in  electric  propulsion  thrusters  has  the  potential  to  significantly  improve  the  design  process  of  these 
devices.  This  goal  has  been  partially  realized  with  numerical  simulations  playing  an  increasing  role  in  the 
design  particularly  of  ion  thrusters  and  Hall  thrusters.  Nevertheless,  the  accurate  simulation  of  many  electric 
propulsion  thrusters  remains  a  significant  challenge  due  to  the  difficulties  associated  with  the  need  to  model 
the  wide  variety  of  physical  phenomena  that  occur  in  these  systems  including  viscous  gas  dynamics, 
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nonequilibrium  collisions,  rarefied  plasma  physics,  and  surface  interactions.  In  terms  of  predicting  thruster 
performance,  accuracy  of  the  existing  numerical  methods  ranges  from  within  5%  for  resisto-jets  and  ion 
thrusters,  to  within  10-  20%  for  arcjets  and  Hall  thrusters. 

There  are  two  main  directions  for  future  work  to  continue  to  improve  the  numerical  modelling  of  electric 
thrusters.  First,  the  numerical  methods  themselves  must  be  improved  in  terms  of  their  physical  accuracy  and 
their  computational  speed.  The  level  of  physical  accuracy  required  in  the  modelling  of  a  particular  thruster 
depends  on  the  nature  of  the  acceleration  mechanism  employed.  To  make  valuable  contributions  to  the  design 
process,  computational  analyses  are  most  useful  when  they  predict  the  correct  trends  in  solution  times  of  a  few 
hours.  This  latter  aspect  is  in  part  addressed  automatically  by  the  continued  increase  in  computational  power 
provided  by  hardware  and  software  developments.  The  second  main  direction  for  improvement  in  the 
simulations  involves  the  more  accurate  determination  of  physical  parameters  that  are  required  by  the 
numerical  formulations.  Examples  of  such  information  include  sputter  yields  for  new  grid  materials  that 
might  be  used  in  gridded  ion  thrusters,  secondary  electron  coefficients  of  wall  materials  used  in  Hall  thrusters, 
electron  mobility  and  other  transport  coefficients,  and  cross  sections  for  new  propellant  species.  These  data 
are  most  likely  to  be  obtained  computationally  due  to  the  development  of  molecular  dynamics  simulation 
methods,  and  the  expense  and  difficulty  of  performing  laboratory  experiments  to  measure  such  data. 
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Table  1.  Properties  of  Various  Electric  Propulsion  Thruster  Types 


Thruster 

Thrust 

Specific  Impulse 

Number  Density 

Flow  Regime 

(mN) 

(sec) 

(m-3) 
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2-100 
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Continuum/Kinetic 

Arcjet 
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Ion  Thruster 
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1500-5000 

-J 
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oo 

Kinetic 

Hall  Thruster 
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Kinetic 

RTO-EN-AVT-162 


PAPER  NBR- 19 


UNCLASSIFIED 


UNCLASSIFIED 


Simulation  of  Electric  Propulsion  Thrusters 


ORGANIZATION 


I _ I _ _ _ , _ _ _ 1 _ I _ l_l _ I _ I _ I _ I _ L 

0.0  0.5  1.0  1.5 

Axial  Distance  (cm) 

Figure  1.  Mach  number  contours  in  a  hydrogen  resisto-jet  (from  Boyd  et  al.  [8]). 
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Figure  2.  Axial  velocity  profiles  for  a  hydrogen  resisto-jet  (from  Boyd  et  al.  [8]). 
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Figure  3.  Specific  impulse  for  a  hydrogen  resisto-jet  (from  Boyd  et  al.  [8]). 
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Figure  4.  Schematic  diagram  illustrating  DSMC  arcjet  analysis  (from  Boyd  [18]). 
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Figure  5.  Axial  velocity  profiles  in  a  hydrogen  arcjet  (from  Boyd  [18]). 


Figure  6.  Radial  velocity  profiles  in  a  hydrogen  arcjet  (from  Boyd  [18]). 
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Figure  7.  Specific  impulse  for  a  hydrogen  arcjet  (from  Boyd  [18]). 
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Figure  9.  Computational  domain  for  ion  optics  simulation  (from  Boyd  and  Crofton  [42]). 


Figure  10.  DSMC-PIC  analysis  of  current  collection  on  the  acceleration  grid  of  the  UK- 10  ion  thruster  (from 

Boyd  and  Crofton  [42]). 
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Figure  11.  Erosion  pattern  around  a  single  aperture  on  the  downstream  face  of  the  acceleration  grid  of  the 
NSTAR  thruster:  photograph  (left)  and  simulation  (right)  (from  Wang  et  al.  [44]). 
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Figure  12.  Erosion  profile  between  the  centers  of  two  adjacent  apertures  for  the  acceleration  grid  of  the 

NSTAR  thruster  (from  Wang  et  al.  [44]). 
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Figure  13.  Schematic  diagram  of  a  Hall  thruster. 
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Figure  14.  One-dimensional  PIC  analysis  of  the  performance  of  an  SPT  Hall  thruster  using  xenon  and 

krypton  propellants  (from  Garrigues  et  al.  [49]). 


RTO-EN-AVT-162 


PAPER  NBR- 27 


UNCLASSIFIED 


UNCLASSIFIED 


Simulation  of  Electric  Propulsion  Thrusters 


ORGANIZATION 


Domain  Exit 


Figure  15.  Computational  domain  for  Hall  thruster  simulation. 


Figure  16.  Centerline  profiles  of  plasma  potential  for  the  P5  Hall  thruster  (from  Koo  and  Boyd  [56]). 
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Figure  17.  Wall  erosion  profiles  for  the  SPT-100  Hall  thruster:  inner  wall  (upper));  outer  wall  (lower)  (from 

Yim  et  al.  [74]). 
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